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MEMORANDUM  REPORT 

SUBJECT:  Models  of  a  Service  System  for  Production  Machine  Maintenance 


1.  Reference: 

References  are  designated  by  bracketed  numbers  and  are  included  in  the  foot¬ 
notes.  The  references  are  also  listed. 

2.  Background 

This  memorandum  reports  on  work  associated  with  a  larger  study  [1],  concerned 
with  the  productivity  of  manufacturing  systems.  A  stochastic  simulation,  called 
TANDEMT,  was  developed  for  the  Ref  [1]  study.  Part  of  TANDEMT  characterizes  a  main¬ 
tenance  system  for  adjusting  and  repairing  machines  on  a  production  line.  To  verify 
this  part  of  the  simulation,  analytic  studies  of  service  systems  were  performed.  The 
models  in  this  study  draw  from  and  extend  classical  queueing  theory,  such  as  in  Ref 
[2].  Whereas  the  analytic  models  were  used  to  verify  a  simulation,  their  use  is  not 
limitec  to  that  purpose.  It  has  been  found  that  one  can  forego  the  simulation  and 
use  only  analytic  models  in  designing  a  machine  maintenance  system*.  Recognizing 


[1]  DD1498,  HQ,  US  ARRCOM,  DRSAR-SA,  March  1983,  title:  Manufacturing  Productivity 
Study. 

[2]  Gross,  Donald  and  Harris,  Carl,  Fundamentals  of  Queueing  Theory,  John  Wiley 
ar-d  Sons,  c.  1974. 

The  analytic  models  are  included  in  a  subroutine  which  can  be  called  from 
TANDEMT.  This  arrangement  is  convenient  since  these  programs  share  input 
data.  In  using  the  routine  for  system  design,  the  minimum  number  of  repair¬ 
men  is  selected  which  satisfies  constraints  on  system  performance  measures. 


the  general*  importance  of  analytic  models  of  this  sort,  I  have  chosen  to  present 
certain  derivations  and  results  of  interest  here. 

3.  Objectives 

One  purpose  of  this  memorandum  is  to  derive  some  results  which  may  be  useful 
in  modeling  a  maintenance  system  which  serves  a  finite  population  of  machines. 
Another  objective  is  to  present  parametric  studies  using  these  models.  A  comparison 
of  results  of  alternative  models  will  also  be  made.  Of  interest  to  analysts  are 
some  numerical  approximations. 

4.  Scope 

In  the  language  of  queueing  theory  the  (finite)  population  of  machines  are 
"customers"  and  the  repairmen  who  perform  maintenance  are  the  "servers". 

Similarly,  the  failures  of  machines  are  regarded  as  demands  for  service,  which 
if  unmet  are  placed  in  a  "service  queue".  I  will  use  these  terms  to  suggest  that 
these  models  are  not  limited  to  a  production  setting.  The  "system"  of  interest 
here  is  the  machine  maintenance  system,  which  is  characterized  by  the  number  of 
repairmen  (or  servers)  and  by  the  probability  distribution  of  service  time. 

These  two  factors  are  the  distinguishing  elements  between  models  in  this  memo¬ 
randum.  Note  that  no  distinction  is  made  between  types  of  repairmen.  For  this 
application  any  repairman  can  render  the  requisite  service.  These  repairmen  are 
not  specialists.  Further,  they  are  assumed  to  work  as  individuals,  rather  than 
in  teams. 


*  Other  applications  of  the  models  considered  here  are  manifold.  Finite 
queueing  models  have  been  used  in  describing  the  process  by  which 
artillery  targets  are  detected  and  "served"  and  in  describing  the  FM 
communications  process  in  tactical  fire  control  nets. 
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5.  Ord-'narily,  the  random  arrival  of  a  "demand"  for  service  is  regarded  as 
occuring  at  a  rate  equal  to  the  product  of  a  unit  rate  times  the  number  of  customers 
outside11'  the  service  system.  The  time-to-next-demand  is  conditionally  an  exponential 
random  variable.  For  most  mature  manufacturing  systems  operating  in  steady-state, 
this  is  a  valid  description  of  demand  for  maintenance.  This  demand  description 

is  used  in  all  models  in  this  report.  Actual  maintenance  experience  often  shows 
that  the  time  to  complete  a  maintenance  action  (time  to  repair)  is  a  lognormal 
random  variable.  To  describe  a  lognormally  distributed  event  in  terms  of  Markov 
processes  involves  approximation.  For  service  times  which  have  a  coefficient  of 
variation  near  unity,  an  exponential  approximation  is  used.  However,  if  the 
standard  deviation  of  service  time  is  significantly  smaller  than  the  mean--as 
is  often  the  case — an  Erlang  distribution  can  be  used  to  describe  the  service 
time.  The  Erlang  distribution,  or  gamma  distribution  with  integer  shape  para¬ 
meter  p  can  be  represented  as  the  distribution  of  a  sum  of  p  random  variables 
each  of  which  is  exponential.  For  example,  in  the  case  of  an  Erlang  distribution 
with  shape  parameter  2  and  mean  x,  an  equivalent  stochastic  model  consists  of  two 
serial  elements  both  of  which  are  exponential  with  common  mean  x/2.  Generally, 
the  coefficient  of  variation  of  a  gamma  distribution  with  shape  parameter  p  is 
1//  p  .  Thus,  if  the  standard  deviation  of  service  time  is  about  0.7  times  the 
mean,  an  Erlang  distribution  with  shape  parameter  2  would  be  an  appropriate 
statistical  model.  The  models  considered  here  have  either  exponential  or  ganma(2) 
service  times. 

6.  Birth-Death  Processes 

Consider  a  service  system  having  a  single  server  whose  service  times  are 
gamma (2)  i.e.,  are  random  variables  from  a  gamma  probability  distribution  with 

shape  parameter  2.  This  is  sometimes  written  r(2).  Let  the  customer  population 
be  m  machines.  Thus,  at  most  m  machines  may  reside  in  the  maintenance  system 
with  at  most  m-1  in  the  queue.  Machine  failures  occur  at  inter-event  times  which 
are  exponential  with  a  rate  (m-k)A,  where  k  machines  are  in  the  system  and  where 


*  The  term  "inside  the  system"  refers  to  residence  in  either  the  service 
queue  or  in  a  service  channel  being  served.  After  being  served  a 
customer  leaves  the  system  and  is  then  "outside  the  system". 
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the  unit  rate  A  is  the  unit  rate  of  entry  or  "birth"  rate  for  a  machine  into 
the  system.  Notationally,  let  MTBF  be  the  mean  time  between  failures  for  the 
population  of  machines.*  Then, 

A  =  1/MTBF.  (1) 

To  model  the  gamma(2)  service  times,  two  stages  of  service  are  assumed,  each 
of  which  is  exponentially  distributed  with  rate  parameter  p.  After  the  second 
stage  of  service  is  completed  the  machine  leaves  the  maintenance  system.  Thus, 
if  Pj  is  the  rate  of  exodus,  or  "death"  rate,  from  the  system 

p  =  2pj  •  (2) 

With  a  mean  time  to  repair  MTTR, 

Pj  =  1/MTTR  .  (3) 

7.  Definitions  of  States  for  a  Single-Server  System 

To  characterize  the  possible  states  of  the  system,  define  and  number  event 
E.  as  follows: 

Ej  =  2j-l  ,  j>l  ,  (4a) 

with  j-1  customers  in  the  queue,  where  the  customer  being  served  is  in  the 
first  stage  of  service;  and 

Ej  =  2j  ,  (4b) 

with  j-1  customers  in  the  queue,  where  the  customer  being  served  is  in  the 
second  stage  of  service.  Eg  denotes  the  null  state.  Notationally,  let 

Pn(t)  =  P{system  is  in  state  En  at  time  t>,  0£  n  _<  2m  .  (5) 

The  object  of  the  following  derivation  is  to  show  how  Pn(t)  is  obtained,  for 
the  admissible  values  of  n,  when  time  becomes  large,  i.e.,  when  the  system  is 
in  stochastic  steady  state. 

8.  State  Transitions  (Single  Server) 

It  is  helpful  to  examine  the  ways  in  which  the  system  can  transition  between 
states.  Figure  1  displays  the  states  as  nodes  and  the  transitions  as  directed 
arcs.  This  state  transition  diagram  also  shows  the  Markov  rate  constant  for 

*  See  Annex  B  for  treatment  of  a  heterogeneous  population  of  machines. 
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Customers  in 
the  System 


Figure  1  .  State  Transition  Diagram  for  a  Single-Server  Service  System  With  Erlang  (2)  Service  Times 


each  indicated  transition  as  a  label  on  the  arc.  Thus,  the  conditional  probability 
of  transition  from  state  0  to  state  1  in  an  incremental  time  interval  h  is  mAh. 
Similarly,  given  the  system  is  in  state  1,  a  transition  to  state  2  in  time 
increment  h  occurs  with  probability  yh,  and  a  transition  to_  state  3  in  time 
increment  h  occurs  with  probability  (m-l)xh.  Following  the  standard  procedure 
for  Markov  processes,  one  writes  the  probability  of  being  in  state  n  at  time 
t+h  in  terms  of  the  probabilities  of  occupying  all  of  the  states  at  time  t. 

For  example. 


P(Eq  at  t+h}  =  P{Eq  at  t}  • 

P{no  transition  in  (t,  t+h)} 

+  PIE^  at  t}P{trans.  from  E^  to  EQ  in  (t,  t+h)}  . 

(6) 

In  our 

notation. 

P0( t+h )  =  pQ(t)  (1  -  mAh)  +  p2(t)yh  . 

(7) 

Then, 

-ft  tP0(t+h)  -  pQ(t)]  =  -mAp0(t)  +  yp2(t)  . 

(8) 

Taking 

the  limit  as  h  approaches  zero  yields 

pQ(t)  =  -mAp0(t)  +  yp2(t)  . 

(9) 

Note  that  this  result  can  also  be  obtained  by  inspection  from  the  state 
transition  diagram.  To  form  the  general  right  hand  side:  All  rates  from  ; 
the  k  th  node  are  summed;  and,  the  negative  of  this  sum  is  the  coefficient 
of  p^.  All  other  terms  of  the  r.h.s.  have  positive  coefficients  corresponding 
to  the  rates  on  arcs  entering  the  k  th  node.  The  j  th  coefficient  multiplies 
p - ( t ) ,  where  j  is  a  node  leading  directly  into  the  k  th  node. 

J 

9.  Kolmogorov  Equations  for  a  One-Server  System 

Following  the  above  procedure  for  all  states  gives  rise  to  the  differential - 
difference  (Kolmogorov)  equations  which  describe  this  system.  Suppressing  the 
notation  for  time  dependence,  one  has: 
pQ  =  -mAp0  +  yp2  . 

Pj  =  mAp0  -  (y+(m-l)A)p1  +  yp4  . 
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For  2  _<  j  <  m-1  , 


p2j-l  “  (m_J+l)A  P2 j_ 3  “  ( ( m- j ) A ) j _ 2  +  uP2j+2 
P2j-  =  (m-j+1 ) A  p2j._2  +  pp2j_1  -  (y+(m-j)A)p2j  . 


And,  for  the  last  two  states, 


p2m-l  Ap2m-3  "  yp2m-l 

p2m  =  Ap2m-2  +  yp2m-l  "  yP2m  ' 


(10) 


In  stochastic  steady  state,  the  time  derivatives  of  all  probabilities  are 
zero.  Setting  dpk/dt  to  zero  for  all  k  in  (10)  yields  a  matrix-vector 
equation  for  the  state  probability  vector  £  ,  with 


The  2m+l  equations  in  the  set  (10)  contain  one  redundant  equation  since  the 
probabilities  are  related  by 


(ID 


One  of  the  equations  can  be  deleted— say,  the  first— and  (11)  substituted  for 

pQ  in  the  second  equation  (the  only  other  place  in  which  pQ  appears).  This 

process  yields 
2m 

-mAEj=l.  p.  -  (y  +  (m-l)A)p1  +  UP4  =  -mA 

(m-l)Ap1  -  (y  +  (m-2) A)p3  +  yp6  =  0 

(m-l)Ap2  +  yp3  -  (y  +  (m-2)A)p4  =  0  , 
and  so  forth. 

The  equivalent  matrix-vector  formulation  is 


A  £  =  J3 


(12) 


with  A  (2m  x  2m)  and  £  and  b  (2m  x  1),  where 
b_'  =  (-mA,  0,  0,  . . . ,  0)  , 


•  •  •  J 


(12b) 


and 


(12c) 
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The  elements  of  A  are: 


aU  “  “mA  “  (y  +  (m-lA) 

alj  =  ~mA  »  2  £  J  ¥  4  £  2m 

a14  =  y  -  mx 
a21  =  y 
a22  =  + 

a2j  =  0  »  J  >  2  • 

For  the  remaining  non-zero  terms 

^aij^  *  2  <.  k  <  m-1  , 

a2k-l,  2k-3  =  (m-k+l) A 
a2k-l,  2k— 1  =  '  (^+(ni-k)A) 

a2k-l,  2k+2  =  u 
a2k,  2k-2  =  (m-k+l)A 

a2k,  2k- 1  =  u 

a2k,  2k  =  “  Cu+(m-k)A)  ; 

a2m-l,  2m- 3  =  x 
a2m-l,  2m- 1  =  -,J 
a2m,  2m- 2  =  A 
a2m,  2m- 1  =  v 

a2m,  2m  =  _,J  •  (13) 


if 
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(14) 


Since  A  is  of  full  rank,  £  may  be  obtained  from 
£  =  A-1b  . 

The  probability  of  the  null  state  is  obtained  from  (11). 

Then,  the  probability  of  the  maintenance  system  having  k  customers  in 
steady  state  is  given  by 

\  =  p2k-l  +  p2k  »  l<k<m.  (15) 


10.  Statistical  Properties  of  the  System 

From  this  probability  distribution  one  obtains  the  mean  and  variance*  of 
the  number  of  customers  in  the  system. 


E[Nsys]  =  Nsys  ”  Ek=l  k7r 
„m  ,  2 


V(Nsys>  =  1  k‘  \  -("sys)2  • 

The  (mathematically)  expected  number  of  customers  in  the  queue  (N^  )  is 
obtained  from 


(16) 

(17) 


E[Nq]  Ek=2 
or 

E[V  -  V  -  <H>0>  5  (18) 

The  variance  of  the  number  in  the  queue 

V[Nq]  =  e”=2  (k-l)\  -  (Nq)2  .  (19) 

Related  properties  are  the  mean  time  a  customer  spends  in  the  system  (W)  and 
the  mean  time  spent  in  the  queue  (W^).  For  all  service  systems  the  value  of 
W  is  greater  than  Wq  by  the  mean  service  time.  In  this  case, 

A  =  Wq  +  V*1!  •  (20) 


*  Operator  notation  is  used  in  displaying  the  mean  and  variance.  The 
expectation  operator  is  denoted  EM,  and  the  variance  operator  V  [  -  ] . 
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Given  W,  equation  (20) is  used  to  obtain  W  . 

With  exponential  arrivals  and  a  finite  customer  population  the  mean  arrival 
rate  is  X(m-N  ) .  Using  this  result  in  Little's  formula  (p.  120,  [2])  yields 

ojr  j 


Nsys 

from  which 


x(m-N  )W  , 
v  sys7 


W  =  N  / [A(m-N  )]  . 
sys  L  v  sys7 


(21) 


Then,  can  be  obtained  from  (20). 
version  of  Little's  formula: 

wq  =  Nq/[*(<n-Nsys))  . 


Alternatively,  is  given  by  another 


(22) 


11 .  Conditional  State  Probability  Distribution 

The  following  derivation  is  to  obtain  the  probability  distributions  of 
waiting  time  in  queue  and  in  the  system,  given  an  arriving  customer  must 
queue.  In  doing  this,  we  first  obtain  the  conditional  probability  distribu¬ 
tion  that  the  system  is  in  the  n  th  state,  given  an  arrival .  Denote,  q^  as 
the  (steady-state)  probability  that  an  arrival  finds  the  system  in  the  n  th 
state.  From  Bayes'  theorem  ,  q  is  proportional  to  the  product  of  two 
probabilities:  (a)  Plthat  an  arrival  occurs,  given  that  the  system  is  in 
state  n}  and  (b)  P{system  is  in  state  n).  Thus, 

q0  -  mx  pQ 
qj_  a  (m-l)x  p1 
q2  «  (m-l)x  p2 
q3  «  (m-2)X  p3 
q4  «  (m-2)X  p4 


q2m.3  -  (m-l)x  p2m_3 
q2m-2  “  P2m-2  ' 

[2]  Gross,  D.  and  Harris  C.  Op.  Cit.  ,  1974. 
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(24) 


Generally,  denoting  the  integer  portion  of  x  by  int[x]  , 
qn  =  (m-int[(n+l)/2])pn/K  ,  0  <  n  <  2m-2  , 

where  the  normalization  constant  K  is  given  by 

K  =  mp0  +  z^‘2  (m  -  int[(n+l)/2])pn  (25) 

The  conditional  probability  density  that  an  arrival  finds  k  other  customers 
in  the  system  is  denoted  by  p^: 

p0  =  q0 

pk  =  q2k-l  +  q2k  ’  l<k<mrl.  (26) 

This  function  is  analogous  to  the  unconditional  probability  density  ir^, 
given  by  equation  (15). 

12.  Notationally,  denote  F  (t)  as  the  distribution  of  waiting  time  in 
queue,  T  : 

q 

Fq(t)  =  P{Tq  <t)  .  (27) 

The  expected  value  of  time  in  queue  from  this  distribution  has  been  denoted 
above  as  W  .  Denote  the  conditional  expectation  of  time  in  queue,  given  a 
queue  as  Wqjw-  For  a  single-server  system,  queueing  occurs  whenever  the 
system  is  not  in  the  null  state.  Therefore, 

Wq|w  =  V(1"q0}  •  (28) 

An  alternative  —  to  equation  (22) — means  of  calculating  W  in  this  case  is 
the  following.  This  expression  is  based  on  the  fact— as  observed  in  the  state 
transition  diagram  (Figure  1)— that  an  odd  state  with  k  customers  (E2k  ,) 
involves  2k  stage  transitions  of  service  to  leave  the  system;  and,  an  even 
state  with  k  customers  (E£k)  involves  2k- 1  stage  transitions  of  service  to 
leave  tie  system.  Each  service-stage  transition  requires  an  expected  time 
of  1/u.  Thus, 

(2kq2k-l+  (2k-l)qzk)  •  (29) 
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Altho  it  is  not  obvious  in  this  case,  the  mean  of  waiting  time  in  queue  given 
by  (29)  equals  the  value  given  by  the  expression  (22)  from  Little's  formula. 
Clearly,  equation  (22)  is  preferable  computationally;  however,  development 
of  (22)  has  motivated  the  method  for  deriving  F  (t). 


13.  Define  an  integer  n*  such  that 
n*  =  n+1  for  n  odd 

n*  =  n-1  for  n  even  , 

with  n  an  index  over  the  system  states:  1  <  n  <  2m-2  . 
Then , 


(30) 


2m  2 

Fq^)  =  qo  +  zn=l  qnPtn*  stages  of  service  are  completed  in  time  <_  t, 

given  an  arrival  finds  the  system  in  the  n  th  state}  .  (31) 

Since  each  stage  of  service  requires  an  exponentially-distributed  time, 
with  mean  u  ,  the  time  to  complete  n*  stages  has  a  cumulative  distribution 
which  is  the  n*-fold  convolution  of  this,  namely  the  gamma  cumulative 
distribution  function  (c.d.f. ): 


Using  (32),  one  can  write  (31)  as 

Fq(t)  =  %  +  ^  •  (33) 

With  a  change  in  the  summation  index  and  using  (30),  equation  (33)  becomes 


Fq(t)  =  %  +  (q2k-l  G(2k,y,t)  +  q2k  G(2k-l,y,t))  .  (34) 

A  computationally  useful  form  of  (34)  may  be  derived  using  Molina's  theorem 
([3]  and  [4],  p.  86): 


rX 

‘o 


un_1  e"u 

(n-1).' 


du  =  1  - 


n-1  xJ  e~X 
j=0  j!  ; 


(35a) 


=  G(n,l,x)  .  (35b) 

[3]  Molina,  E.C.  Poisson's  Exponential  Binomial  Limit,  Van  Nostrand,  c.  1942. 

[4]  Schlenker,  G.  Numerical  Methods  in  Renewal  Theory,  (AD  828276),  Hdqts 
USA  WECOM,  February  1968. 
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Molina's  theorem  relates  a  standardized  gamma  c.d.f.  to  a  sum  of  Poisson  terms. 
Note  that  the  c.d.f.  in  (32)  with  parameters  n*  and  y  and  with  argument  t  is, 
via  a  change  of  variable,  equivalent  to  G(n*,l,yt).  Thus,  the  gamma  c.d.f.'s 
in  (34)  can  be  replaced  by  the  sums  required  by  (35)  using  an  argument  x, 
where 


x  =  yt  . 

This  substitution  yieldsthe  formula  for  calculating  F  (t) 

J  --x  q 

(r,  / 

'k=l 


F  (t)  =  1  -  zm_1  (a  z2^1  —  e"A  +  n  TZk~2  xJ  e_Xl 
1  H=i  ^-:=n  j|  +  q2k  ^j=Q  J1  ' 


’2k-l  j=0 


(36) 

(37) 


or 


Fq(t)  =  1  -  q1  (e  x  +  xe"x)  -  q2e 


-x 


-  e 


-_x  ym~l  In  v2k~l  xJ  +  n  v2k-2  xJ  \ 

Zk=2  q2k-l  Ej=0  j!  +  q2k  Ej=0  jl  ]  * 


(38) 


This  unconditional  probability  distribution  for  time  in  queue  can  also  be 
used  to  get  the  conditional  c.d.f.,  given  an  arrival  must  queue.  The  latter 
is  denoted  by  Fq|w(t).  Since  qQ  is  the  probability  that  an  arrival  does  not 
need  to  queue, 


Fq|w(t)  =  (Fq(t)  "  q0)/(1_q0)  *  <39> 

Using  (34)  and  (39), 

Fq|w(t)  =  1^  Ek=l  (q2k-l  +  q2k  G(2k-l,y,t))  .  (40) 

Note  that  the  conditional  c.d.f.  of  time  in  queue  is  a  weighted  sum  of 
gamma  distributions. 


14.  Mean  and  Variance  of  Conditional  Waiting  Time  in  Queue 

The  j  th  origin  moment  of  G(n,y,t)  is  given,  from  p.  33,  [4],  as 

aj  =  y"J  n^=1  (n+i-1)  . 


(41) 


[4]  Schlenker,  G.  Op.  Cit. 
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Specifically,  the  first  and  second  origin  moments  are: 

.  _  n  -1 
aj  -  ny 

a2  -  n(n+l)y  2  .  (42) 

Using  these  results,  one  can  obtain  expressions  for  the  first  and  second 
origin  moments  for  the  conditional  time  in  queue  (T  ,  ). 

From  (40),  q,W 

E[Tq|w'  "  l-qQ  Ek=l  q2k-l  2kp  +  (2k'1)p  1  ■  (43) 

and 

E^q|w^  1-q^  Ek=l  ^2k-l  ^k(2k+l)y  +  q^  2k(2k-l)y  ^  .  (44) 

Then,  the  variance  of  the  conditional  time  in  queue  is  obtained  from 

V|-Tq|w]  “  E|-Tq|w^  "  ^E|-Tq|w^  *  (45) 

The  expression  E[Tqjw]  was  earlier  denoted  by  W^,  and  was  related  to  W 
by  (28).  Thus,  equation  (28)  together  with  (22)  provide  a  consistency  check 
for  the  value  given  by  (43).  This  check  is  performed  in  the  attached  computer 
program  FINITE. ME2.Q. 

15.  Multiple  Servers  With  Exponential  Service  Times 

The  formulas  derived  at  this  point  pertain  to  the  case  in  which  a  single 
server  exhibits  service  times  which  are  Erlang(2).  The  Kolmogorov  equations 
for  this  case,  tho  simple  enough,  did  not  lead  to  a  solution  of  the  state 
probabilities  (pn)  in  closed  form.  All  equations  involved  numerical  solution 
with  a  degree  of  complexity  requiring  computer  assistance.  Another  model  of 
a  finite-customer  population  which  does  lend  itself  to  a  somewhat  simpler 
solution  is  one  to  which  I  now  turn.  In  this  model  the  service  times  are 
exponential  random  variables— i  .e. ,  Erlang(l)  distributed.  Further,  the 
number  of  servers  is  generalized  to  an  arbitrary  integer  c.  This  model  has 
been  developed  in  some  detail  by  Gross  and  Harris  [2].  For  this  model  p 

n 

represents  the  probability  that  n  customers  are  found  in  the  service  system 
in  stochastic  steady  state.  From  p.  118  of  Ref  [2], 

[2]  Gross,  D.  and  Harris,  C.  Op.  Cit. 
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0  <  n  <  c 


(46a) 


-  /HU 
qn 


n! 

,n-c  ,  y  P0 


c  <  n  <  m 


(46b) 


The  value  of  the  probability  for  the  null  state  is  found  by  setting  the 
sum  of  the  state  probabilities  to  unity. 


zc-l  (mj  (X}n  +  m  (m} 
n=0  n  vy;  n=c  vrr 


r\l 


n-c  , 

c  c. 


<7>" 


(47) 


Equations  for  the  mean  and  variance  of  the  number  of  customers  is  the 
system— analogous  to  (16)  and  (17)--are: 

V  =  Ek=l  k  pk  (48) 

V<Nsys>  ■  k'  »k  '  <%s)2-  (49) 

Since  queueing  occurs  when  N  >c,  the  expected  number  of  customers  in 
queue  is  given  by 

"q  =  (k-c)  pk 

or 

"q  =  Ssys  ‘  c  +  (c-k>  pk  •  <5°) 

Little's  formula  is  also  valid  for  this  model;  hence,  the  mean  system 
waiting  time  and  mean  waiting  time  in  queue  are  given  by  (21)  and  (22), 
respectively. 


16.  Distribution  of  Time  in  Queue  (Exponential  Services) 

The  notational  conventions  used  in  pgf  11  for  the  model  with  Erlang(2) 
service  times  are  also  used  here.  The  conditional  probability  that  an 
arriving  customer  finds  the  system  in  the  n  th  state  is  denoted  q  .  Bayes1 
rule  recuires  that 

qn  =  (m-n)Pn/K  »  (51) 

Where  the  normalization  constant  K  is  found  by  requiring  that 

jn-1  „  _  . 

En=0  qn  “  *• 

Thus, 

K  =  (m-n)pn  .  (52) 
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The  following  derivation  for  the  c.d.f.  of  time  in  queue  for  a  c-server 
system  is  a  generalization  of  a  result  found  on  p.  125  of  [2]. 

Let  F  (t)  be  the  cumulative  (unconditional)  distribution  of  waiting  time 

M 

in  queue  as  in  (27).  Then, 


Fq(t)  =  qnP{n-c+l  job  completions  occur  in  time  <_  t,  given 
an  arrival  finds  n  in  the  system}  +  P{system  is  open)  ,  (53a) 


where 


c  —  i 

P{system  is  open}  =  En=Q  qn  .  (53b) 

In  this  model  a  job  completion  occurs  in  exponentially  distributed  time 
with  rate  constant  cy,  when  the  system  is  full. 

Therefore,  the 

P{n*  completions  occur  in  time  <_  t} 
is  given  by 

G(n*,  cy ,t)  , 

with  the  gamma  c.d.f.  defined  by  (32). 

Thus, 


Zn=c  qnG^n  c+1’c^’t)  +  p=0  qn  ' 


(54) 


Using  Molina's  theorem  (35),  a  computationally  useful  form  of  (54)  is  obtained. 
F  ( t )  =  1  -  q  e'Cwt  -e‘c,,t  l”’1.,  q  eI"' 

q  Mc  n=c+l  Mn  j=0  j!  .  (55) 

From  (54)  it  is  seen  that  the  unconditional  expected  time  in  queue  is  the 
weighted  sum  of  the  means  of  gamma  distributions  having  different  shape 
parameters.  With  (42), 


E[Tq]  =  qn(n-c«)/<c) 


(56) 


Similarly, 

E[Tq]  =  eJJ:J  qn(n-c+l)(n-c+2)/(cy)2  . 

Then,  the  unconditional  variance  of  time  in  queue  is 
V(Tq)  =  E[T2]  -  (E[Tq])2  . 


(57) 

(58) 


[2]  Gross,  D.  and  Harris,  C.  Op.  Cit. 
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The  conditional  probability  distribution  of  queueing  time,  given  an  arrival 
must  wait,  is  obtained  from  F  (t)  by 

Fq|w(t)  =  (Fq(t)  -  P{open})/(l-P{open})  ,  (59) 

where  P{open}  is  given  by  (53b). 

This  expression  is  analogous  to  (39)  for  the  single-server  model. 

The  first  two  origin  moments  of  the  conditional  probability  distribution 
of  queueing  time  are  obtained  from  the  corresponding  unconditional  moments 
(56,  57)  by  dividing  each  by  l-P(open).  The  variance  of  this  conditional 
c.d.f.  is  obtained  from  the  general  expression  in  (45). 

17.  Numerical  Results 

Before  proceeding  with  the  derivation  of  other  analytic  relationships, 
it  is  interesting  to  apply  the  formulas  developed  at  this  point  to  specific 
examples.  The  numerical  values  chosen  are  motivated  by  the  author's  experience 
with  maintenance  of  machines  for  producing  metal  parts  of  conventional 
ammunition.  In  performing  the  calculations  two  computer  programs  have  been 
useful:  FINITE. ME2.Q  and  BEST. SERVICE.  These  subprograms  are  listed  in 
Annex  A  together  with  an  executive  driver.  The  first  program  implements 
the  model  with  conditional  Poisson*  arrivals  and  Erlang(2)  services.  The 
second  program  calculates  all  the  statistics  for  the  multiserver  model  having 
conditionally  Poisson  arrivals  and  exponential  services.  Both  programs  treat 
the  population  of  customers  as  finite.  In  order  to  contrast  results  from  a 
finite-population  model  with  a  model  having  an  infinite  customer  population, 

I  also  evaluate  the  latter  model  with  an  arrival  rate  parameter  (x^)  chosen 
to  yield  the  maximum  arrival  rate  in  the  finite  models.  Thus, 

Xro  =  mx  .  (60) 


*  Poisson  arrivals  implies  that  the  time  between  arriving  customers  is  an 
exponentially  distributed  random  variable  and  the  arrival  rate  is  a 
constant  parameter.  In  this  case,  I  mean  by  "conditionally  Poisson" 
that  the  arrival  rate  per  customer  outside  the  system  is  a  constant. 


17 


The  model  with  infinite  customer  population  which  is  chosen  for  comparison  with 
the  models  of  this  report  has  an  Erlang(2)  service  time  distribution  with  a 
single  server.  Certain  mean-value  statistics  for  this  case  are  presented  in 
[2],  starting  on  page  163.  This  model  is  referred  to  as  (M/E2/l)  by  Gross 
and  Harris  [2]. 

The  expected  waiting  time--Wq,  in  the  notation  of  [2]--for  (M/E2/l)  is 

3Aoo 

Wq  =  4m  (m  -X  ) 

*  00  oo  oo 


where  ^  is  the  stage  service  rate,  denoted  elsewhere  as  y . 

The  expected  waiting  time  in  the  system  (W)  is  given  by  (20),  since  this  equation 
is  independent  of  population  size.  Little's  formula  yields  the  expected  number 
in  the  queue  (L  )  and  the  expected  number  in  the  system  (L): 

H 


L  =  A  Wn 
q  *  q 

L  =  A  W  . 


(62a) 

(62b) 


The  probability  that  the  service  system  is  empty  (pg)  is  also  given  in  [2]: 

P{N=0}  =  pQ  =  1  -  XjUco  .  (63) 

Equations  (60)  thru  (63)  were  evaluated  for  comparison  with  their  counterparts 
of  the  finite-population  models.  A  numerical  comparison  of  the  infinite- 
population  model  with  the  finite  population  models  is  made  in  Table  1.  The 
notational  designations  for  the  three  models  used  in  Table  1  are  Infinite 
r(2).  Finite  r(2),  and  Finite  Exponential  to  indicate  the  customer  population 
size  and  the  service-time  distribution  function,  respectively.  The  time 
between  arrivals  (in  minutes)  and  the  number  of  customers  in  the  population 
are  treated  as  parameters.  One  notes  that  certain  output  variables  agree 
fairly  well  across  models,  such  as  Pg,  whereas  others,  such  as  W^,  do  not. 

Due  to  the  fact  that  exceeds  the  average  arrival  rate  in  the  finite  models, 
the  average  number  in  the  system  and  in  the  queue  exceeds  their  finite- 
population  counterparts.  However,  it  is  noteworthy  that  when  the  time  between 
arrivals  is  large  (or  X  is  small)  and  the  customer  population  is  large,  the 
values  of  E[N]  and  E[Q]  from  the  infinite  model  approximate  the  corresponding 
values  in  the  finite  r(2)  model. 


18 


TABLE  1 


COMPARISON  OF  THE  STOCHASTIC  PROPERTIES  OF 
SEVERAL  MODELS  OF  MAINTENANCE  SERVICE  SYSTEMS 


Assumptions : 


Distribution  of  time  between  arrivals  is  exponential  for 
each  customer  in  a  finite  population  of  customers. 

There  is  a  single  server  with  mean  service  time  of  20  minutes. 

l/X  (1> 

Number 

(2) 

System 

Type  of  Model  ^ 

(min) 

Customers 

Property 

Inf.,  r (2) 

Finite,  T ( 2 ) 

Finite,  Exp 

120 

2 

E[N] 

0.4583 

0.3128 

0.3200 

SD[N] 

0.5274 

0.5455 

E[Q] 

0.1250 

0.0316 

0.0400 

P{N=0} 

0.6667 

0.7188 

0.7200 

Wq 

7.500 

2.249 

2.857 

Wq  |  W 

14.800 

20.000 

3 

E[N] 

0.8750 

0.5191 

0.5410 

sd[n] 

0.6934 

0.7372 

E[Q] 

0.3750 

0.1057 

0.1311 

P{N=0} 

0.5000 

0.5865 

0.5902 

Wq 

15.000 

5.111 

6.400 

Wq  |  W 

17.181 

22.857 

5 

E[N] 

3.9582 

1.0936 

1.1624 

sd[n] 

1.0479 

1.1507 

E[Q] 

3.1249 

0.4425 

0.5228 

P{N=0} 

0.1667 

0.3489 

0.3604 

Wq 

74.998 

13.593 

16.348 

Wq  |  W 

24.164 

30.820 

10 

E[N] 

OO 

4.1705 

4.2588 

SD[N] 

1.9462 

2.1538 

E[Q] 

3.1989 

3.3020 

P{N=0} 

0.0284 

0.0431 

Wq 

65.849 

69.017 

Wq  |  W 

68.824 

74.625 

20 

E[N] 

00 

14.0000 

14.0000 

SD[N] 

2.1354 

2.4494 

E[Q] 

13.0000 

13.0000 

P{N=0} 

0.0000 

0.0000 

Wq 

260.000 

260.001 

Wq  |  W 

260.000 

260.004 
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TABLE  1  (CONT) 


COMPARISON  OF  THE  STOCHASTIC  PROPERTIES  OF 
SEVERAL  MODELS  OF  MAINTENANCE  SERVICE  SYSTEMS 


Assumptions : 

Distribution  of  time  between  arrivals  is  exponential  for 
each  customer  in  a  finite  population  of  customers. 

There  is  a  single  server  with  mean  service  time  of  20  minutes. 


1A  (1) 

(min) 

Number 

Customers 

(2) 

System 

Property 

Type  of  Model  ^ 

Inf.,  T(2) 

Finite,  F(2) 

Finite,  Exp 

1200 

20 

E[N] 

0.4583 

0.4336 

0.4668 

SD[N] 

0.7210 

0.8012 

E[Q] 

0.1250 

0.1075 

0.1413 

P{N=0} 

0.6667 

0.6739 

0.6744 

Wq 

7.500 

6.591 

8.679 

Wq  |  W 

21.138 

28.049 

30 

E[N] 

0.8750 

0.7998 

0.8928 

SD[N] 

1.0634 

1.2358 

E[Q] 

0.3750 

0.3131 

0.4076 

P{N=0} 

0.5000 

0.5133 

0.5149 

Wq 

15.000 

12.868 

16.806 

Notes  for  Table  1. 


(1)  The  mean  time  between  arrivals  per  customer  is  1/A.  Thus,  the  maximum 
arrival  rate  is  A  times  the  number  of  customers. 

(2)  The  following  notation  is  used  to  denote  system  properties: 


E[N] 

Statistically  expected  number  of  customers  in  the  service 
system. 

sd[n] 

Standard  deviation  of  the  number  of  customers  in  the 
service  system. 

E[Q] 

Expected  value  of  the  number  in  the  service  queue. 

P{N=0} 

Probability  that  the  system  is  empty. 

Wq 

Expected  value  of  the  waiting  time  (minutes)  in  the 
service  queue. 

Wq  |W 

Expected  waiting  time  (minutes)  given  a  customer 
must  queue. 

(3)  Models  are  characterized  here  by  the  population  of  customers  and  by  the 
form  of  the  service-time  distribution.  Thus,  MInf.,  r(2)M  signifies  an 
infinite  population  of  customers  and  a  service-time  distribution  which 
is  gamma  (or  Erlang)  with  shape  parameter  2.  The  form  of  the  c.d.f.  of 
time  t  is 

F(t)  -  1  -  (l+gt)i6* 

with  parameter  g.  The  mean  of  t  with  this  distribution  is  2/g. 

The  other  model  has  a  finite  customer  population  and  exponentially 
distributed  service  times. 
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18.  Comparisons  Between  Finite  Models 

The  model  with  exponential  service  times  is  computationally  much  easier 
than  the  r(2)  model.  Therefore,  it  is  useful  to  ask  for  what  statistics  the 
models  are  in  reasonable  agreement.  For  these,  the  model  with  exponential 
service  will  do.  Some  observations  can  be  made  regarding  the  similarity  of 
results  from  the  finite  r(2)  model  versus  the  finite  exponential  model. 

Generally,  very  good  agreement  exists  between  these  models  for  the  values 
of  E[N],  SD[N],  and  P{N=0}.  Poorer  agreement  exists  for  E[Q],  and  very 
poor  agreement  exists  for  and  W^,w  .  Clearly,  the  difference  in  dispersion 
of  service  times  has  the  greatest  effect  upon  mean  waiting  time  (or  conditional 
mean  waiting  time).  Differences  in  the  conditional  waiting  time  distributions 
are  displayed  explicitly  in  Figures  2  and  3.  These  figures  differ  in  only  one 
parameter-customer  population  size.  Figures  4  and  5  illustrate  the  parametric 
effect  of  population  size  for  the  model  with  exponential  services.  Using  this 
model,  the  parametric  effect  of  number  of  servers  on  the  c.d.f.  of  waiting  time 
in  queue  is  shown  in  Figure  6.  These  probability  distributions  are  plotted  on 
Wei  bull  probability  paper  for  convenience,  since  extreme  variation  in  waiting 
time  is  shown  for  any  value  of  a  parameter.  The  following  interesting  observa¬ 
tion  can  be  made  from  Figure  6.  When  the  number  of  servers  increases— indicating 
improving  service— the  conditional  distribution  function  of  time  in  queue  closely 
approximates  an  exponential  distribution.  An  exponential  distribution  plots  as  a 
straight  line  on  Weibull  paper  with  a  slope  of  unity  (Note  scale  difference  of 
abscissa  and  ordinate.)  Even  tho  the  mathematical  form  of  (55)  is  clearly  not 
that  of  an  exponential  c.d.f.,  nevertheless  an  exponential  form  is  a  good 
approximation  under  certain  conditions.  After  examining  a  variety  of  parametric 
variations  with  this  model,  I  conclude  that  a  condition  under  which  the 
exponential  approximation  is  pragmatically  adequate  is  the  following:  when 
the  mean  value  of  conditional  waiting  time  in  queue  does  not  exceed  the  mean 
service  time.  Examples  illustrating  this  point  are  given  in  Table  2.  The 
exponential  distribution,  with  single  parameter  0,  which  approximates  the 
conditional  c.d.f.  of  waiting  time  in  queue  is  denoted  by  F(t,9).  The  exact 
c.d.f.  is  denoted  F(t). 

F(t,0)  =  1  -  exp(-t/0)  .  (64) 
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Customer  population  =  10. 

Mean  service  time  =  20  minutes. 


The  parameter  0  is  assigned  the  exact  conditional  mean  waiting  time  (W  ). 

This  value  can  be  obtained  from  easily  computed  quantities  as  follows. ^ 

The  system  mean  waiting  time  is  calculated  from  (21).  Then,  w  is  calculated 
via  (20).  Both  of  these  relations  exploit  Little's  formula.  Then,  Bayes1 
theorem  requires 

Wq|w  =  Wq/(1“Pq)  *  (65) 

where  pq  is  the  probability  that  an  arrival  must  queue.  With  c  servers, 

Pq  =  1  "  En=0  qn  *  (66) 

For  a  single  server, 

q0  =  1  -  Pq  (67) 

q0  =  mq0/En=l  (m-n)  P„  '  (68) 

Four  pairs  of  c.d.f.  s  are  shown  in  Table  2.  The  parameter  which  distinguishes 
each  pair  is  the  number  of  servers,  which  varies  from  1  to  4.  The  constants 
across  comparisons  are  customer  population  size,  mean  interarrival  time  per 
customer,  and  mean  service  time.  The  mean  service  time  is  20  minutes. 

Note  that  the  conditional  mean  waiting  time  is  less  than  20  minutes  in 
the  last  three  comparisons,  but  not  the  first.  In  the  first  case,  with  a 
single  server  the  approximation  is  poor;  whereas,  in  the  other  cases  the 
approximation  is  quite  good. 

19.  Concl us  ions 

The  inferences  drawn  from  the  numerical  examples  are  summarized  here. 

o  Unless  the  activity  parameter  (mA/y)  is  smaller  than  about  1/3, 
use  of  queueing  theory  for  an  infinite  customer  population  to 
describe  the  behavior  of  a  finite  queueing  model  is  quite  inexact. 
Further,  there  are  no  compelling  numerical  reasons  for  using 
results  applicable  to  infinite  systems. 

o  Numerical  problems  in  calculating  the  system  steady-state  probability 
vector  were  non-existent  using  double-precision  arithmetic  on  the 
-’RIME  750  minicomputer. 

o  If  expected  values  of  service  system  properties  such  as  number  of 
customers  in  the  system  and/or  in  the  queue  are  of  primary  concern, 
the  differences  between  finite  models  having  exponential  and  r(2) 
service  times  seem  insignificant  in  practice. 
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o  One  must  qualify  the  previous  conclusion.  The  form  of  the  probability 
distribution  of  service  time  should  be  selected  carefully  when  the 
conditional  mean  time  in  queue  must  be  accurate.  Even  with  the  same 
mean,  different  service  time  c.d.f. 's  may  yield  quite  different  mean 
waiting  times.  This  observation  holds  for  exponential  versus  r(2) 
service-time  distributions. 

o  The  functional  form  of  the  conditional  probability  distribution  of 
waiting  time  in  queue  is  generally  quite  complex.  But,  a  simple 
approximation  exists  under  certain  conditions.  For  a  model  with 
exponential  service  times,  the  approximation  is  valid  when  the  conditional 
average  waiting  time  in  queue  does  not  exceed  the  mean  service  time. 

20.  Recap  and  Survey  of  Other  Derivations 

A  review  at  this  point  is  appropriate.  I  have  discussed  two  mathematical 
models  used  to  describe  the  maintenance  system  for  production  machines.  Because 
the  population  of  machines  is  finite,  these  queueing  models  treat  a  finite 
customer  population.  In  both  cases  machine  failure  is  assumed  to  occur 
stochastically  at  a  constant  rate  during  machine  operation.  A  failure 
constitutes  the  birth  of  a  service  demand.  Similarly,  the  completion  of  service 
constitutes  the  death  of  the  service  action.  With  this  perspective  a  Markov 
birth-death  process  was  used  to  describe  a  particular  service  system  probabil isti - 
cally.  This  single-server  system  completes  service  in  a  time  that  has  an 
Erlang(2)--  or  r(2)--probabil  ity  distribution.  For  contrast,  a  second  model 
was  described  which  exhibits  exponentially  distributed  service  times.  The 
mathematical  development  of  the  latter  model  was  nearly  complete  in  the 
literature  on  queueing.  Only  the  conditional  c.d.f.  of  waiting  time  in  queue 
needed  to  be  derived  for  the  multiple-server  case.  This  was  done  next.  Without 
further  mathematical  developments,  I  presented  some  numerical  results.  Inferences 
were  made  from  these  results.  In  the  remainder  of  this  memo  I  will  generalize  the 
analytic  results  of  the  finite  r(2)  model  to  treat  several  servers.  This  generaliza¬ 
tion  proceeds  in  steps.  First,  two  servers  are  considered.  The  system  state- 
transition  diagram  is  constructed.  From  this  the  Kolmogorov  equations  are  written. 
Then  the  method  of  numerical  solution  is  indicated  and  certain  statistical 
properties  are  derived.  In  deriving  these  results,  equations  for  the  general 
multiple  server  case  are  presented.  Finally,  a  finite  r(2)  model  with  three 
servers  is  considered.  The  same  procedure  for  deriving  results  is  followed  here 
as  was  used  for  the  one-  and  two-server  cases. 
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21.  Two  Servers  With  Erlanq(2)  Service 

The  equations  derived  in  paragraphs  6  thru  14  above  apply  to  a  service 
system  with  a  single  server.  This  case  was  selected  initially  because  of 
its  mathematical  simplicity.  Frequently,  several  repairmen  work  independently 
and  in  parallel  on  machines  of  (a  segment  of)  a  production  system.  To  introduce 
the  treatment  of  multiple  servers,  consider  the  case  where  2  servers  are  present. 
Each  is  assumed  to  have  an  identical  service  time  distribution  which  (again)  is 
Erlang(2).  The  derivation  of  results  follows  the  pattern  used  for  a  single 
server,  with  appropriate  modification  in  the  definition  of  system  states.  In 
labeling  the  states  of  the  system,  it  is  convenient  to  start  with  a  notation 
which  uses  double  subscripts.  Later,  the  doubly  subscripted  state  variable  is 
replaced  with  a  variable  having  a  single  subscript.  The  procedure  followed  in 
this  case  is  easily  generalized  to  cases  having  more  servers.  An  identical 
approach  was  used  in  an  analytic  communications  model  for  artillery  with  CLGP 
Copperhead  (p.  7,  [5]). 

22.  I  he  state  of  a  two-server  system  is  characterized  by  the  number  of  customers 
in  the  system  and  by  the  stage  of  service  of  each  of  the  servers.  The  first 
server  is  arbitrarily  considered  to  take  the  first  customer.  Thus,  when  only 
one  customer  is  in  the  system,  the  first  server  is  either  in  the  first  stage 

of  service  or  the  second,  with  the  other  (2nd)  server  idle.  This  situation  can 
be  represented  iconically  by  the  two  system  configurations: 

3  10  or  0  0  1 

0  0  0  0  , 

where  the  first  integer  of  the  first  row  is  the  number  of  customers  in  queue 
and  where  a  1  is  placed  in  the  active  stage  of  service  for  each  server. 

States  a~e  symbolized  by  the  label  E. . ,  where  the  index  i  equals  the  number 
of  customers  (k)  in  the  system  and  the  second  index  j  takes  on  values  from 
1  to  3  indicating  service  configuration.  All  of  the  state  configurations 
and  indices  are  shown  in  Table  3. 


[5]  Sch^enker,  G.  "Field  Artillery  Communication  Studies  with  Application 

to  CLGP,"  Systems  Analysis  Directorate  Activities  Summary  for  April  1977. 
May  1977. 
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TABLE  3 


STATES  OF  THE  SERVICE  SYSTEM  WITH 
TWO  SERVERS  AND  TWO  STAGES  OF  SERVICE 


Label 

(i  >j) 

Q 

Configuration 

Stage  1  Stage  2 

Customers  in 

the  System 

0 

0 

0 

0 

0 

0 

0 

1.1 

0 

1 

0 

1 

0 

0 

1,2 

0 

0 

1 

1 

0 

0 

2,1 

0 

1 

0 

2 

1 

0 

2,2 

0 

0 

1 

2 

1 

0 

2,3 

0 

0 

1 

2 

0 

1 

3,1 

1 

1 

0 

3 

1 

0 

3,2 

1 

0 

1 

3 

1 

0 

3,3 

1 

0 

1 

3 

0 

1 

Generally,  for 

2<_k_<  m. 

k,  1 

k-2 

1 

0 

k 

1 

0 

k,2 

k-2 

0 

1 

k 

1 

0 

k,3 

k-2 

0 

1 

k 

0 

1 
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Customers  in 
the  System 
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i 
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Figure  7.  State  Transition  Diagram  for  a  Two-Server  System  with  Erlang  (2)  Service  Times 


23.  Kolmogorov  Equations  for  a  Two-Server  System 

As  was  done  for  the  single-server  system,  one  can  construct  a  state 

transition  diagram  for  this  system  preparatory  to  writing  the  Kolmogorov 

equations.  However,  the  states  in  this  diagram  are  labeled  with  the 

double-index  notation.  This  diagram  is  shown  in  Figure  7.  We  use  the 

same  notation  for  the  unit  birth  rate  (a)  into  the  system  and  the  unit 

death  rate  (y)  out  of  the  system.  The  definitions  of  A  and  y  are  given 

in  equations  (1)  and  (3).  As  assumed  earlier,  the  time  between  arrivals 

for  customers  outside  the  system  is  exponential.  Since  each  server  is 

assumed  to  act  independently,  the  transition  rate  from  the  E,  state  to 

1)1 

the  state  is  y  (the  stagewise  service  rate  for  one  server).  If.  both 
servers  worked  together,  the  mean  stage  service  time  at  this  point  would 
be  less  than  vi  . 

24.  For  this  system  define 

P-j-jU)  =  P{system  is  in  state  E.  .  at  time  t}  .  (69) 

1  J  1  J 

Then,  using  the  state  transition  diagram,  one  can  write  the  Kolmogorov 
equations  by  inspection.  Notational  dependence  on  time  is  suppressed. 

For  the  non-queueing  states, 


(70) 


For  i  customers  in  the  system  with 
3  <:  i  £  m-1  , 


(71) 
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The  last  three  states  are  described  by 


xP 


m-1,1 


XP 


m 


_1  O  +  2yP  .  -  2yP  9 
m  3 1  m  ^ 


X P  .  o  +  yP  9 
m-1,3  m,2 


(72) 


25.  Steady-State  Equations  (Two  Servers) 

To  obtain  the  steady-state  solution  to  the  above  equations,  all 
derivatives  are  set  to  zero.  At  this  point,  the  double-subscript  notation 
is  changed  to  a  single  subscript  notation  in  denoting  states.  Let  the 
probability  of  the  k  th  state,  in  stochastic  steady  state,  be  denoted  by 
Pk,  so  that 

p0  -  p0<») 

pi  ■  pn(") 

P2  =  P12^°°^ 

p3  =  P21^ 

P4  =  p22(a>) 

P5  ”  p 23^”)  (73a) 

In  general ,  for  2  £  i  £  m  and  1  <  j  <  3  , 

k  =  3 ( i - 1 )  +  j-1 
and 

p3i+:-4  =  Pij(“>  '  73b) 

Note  that  max(k)  =  3m-l. 

With  this  notational  change,  the  steady-state  system  equations  are, 
from  (70,  71,  72,  73)  , 

-mxpQ  +  yp2  =  0 

mXpQ  -  ((m-l)x  +  y)p2  +  yp4  =  0 


35 


(74) 


yPj  "  ( (m-l ) A  +  y)pg  +  2yPj-  =  0 
(m-l)Ap1  -  ( (rn-2 ) A  +  2y)p3  +  yp7  =  0 
(m-l)xp2  +  2yp3  -  ( (rn-2)A  +  2y)p4  +  2ypg  =  0 
yP4  -  ( (m-2 ) A  +  2y)p5  =  0  . 


Since  the  probabilities  {pk)  must  sum  to  unity,  the  first  (redundant) 
equation  in  (74)  can  be  discarded  and  pg  eliminated  from  the  second  via 
the  substitution 


(75) 


This  operation  permits  the  equations  for  {p^}  to  be  written  in  matrix- 
vector  form: 


A  £  =  b^  , 

where 


(76) 


£  and  b^  are  ( 3m- 1  x  1)  and  A  is  (3m-l  x  3m-l) . 
From  the  second  equation  in  (74),  with  (75), 

a^  =  -2mA  +  A  -  y 


a 24  y  —  mA 

alj  =  ~mA  5  1  <  j  <  3m-l  and  j  t  4  .  (77) 

A1  so 

b. '  =  ("mA  >  0,  0,  ...,  0)  .  (78) 

The  elements  of  A~{arc> —  can  be  assigned  in  groups  of  3  rows  each  as 
fol lows: 

For  6  _<  k  _<  3m-4,  corresponding  to  3  _<  i  _<  m-l  and  1  £  j  £  3  in  the  double¬ 
subscript  notation,. 

let  '  r  =  3( i-1)  .  (79a) 

Then, 

ar,3(i-2)  =  ^m_1  +  ^A 
ar,3(i-l)  =  +  2tJ) 

ar,3(i+l)  =  y  .  (79b) 
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Now  assign 


r  =  3(1-1)  +  1  (80a) 

in 

ar,3(i-2)  +  1  =  (m-i+1)A 
ar,3(i-l)  =  2y 

ar  ,3(i-l)  +  1  =  +  2y) 

ar,3i+2  =  2y  •  (80b) 

Finally,  for 

r  =  3( i - 1 )  +  2  ,  (8i) 

ar,3(i-2)  +  2  =  (m-i+1)A 
ar,3(i-l)  +  1  =  y 

ar,3(i-l)  +  2  =  "((m_1>A  +  2p)  •  (82) 

For  the  last  three  rows  in  A, 

a3m-3,3m-6  =  A 
a 3m- 3, 3m- 3  =  "2y 
a3m-2,3m-5  =  A 
a3m-2,3m-3  =  2y 
a3m-2,3m-2  =  _2y 
a3m-l ,3m-4  =  A 
a3m-l,3m-2  =  y 

a3m-l,3m-l  =  _2y  *  (83) 

As  in  the  case  with  a  single  server,  the  state  probability  vector  (q)  is 
obtained  via  (14),  with  ^  and  A  given  in  (78)  thru  (83). 
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26.  Statistics  for  the  Two-Server  System 

To  calculate  the  mean  and  variance  of  number  of  customers  in  the  system 
requires  aggregation  of  the  system-state  probabilities.  Recall  that 
is  defined  as  the  probability  that  k  customers  are  in  the  system.  Referring 
to  Table  3  for  the  definition  of  states, 


*1 


=  P  +  P 
11  *12 


=  zj=i  pkj 


2  <  k  <  m  . 


(84a) 


Using  the  single-subscript  notation,  from  (72  and  73), 
*0  =  p0 


*1  =  Pi  +  P2 

"k  =  p3k-3  +  p3k-2  +  p3k-l  ’  2  -  k  -  m  *  (84b) 

The  mean  and  variance  of  number  of  customers  in  the  system  are  obtained 
from  (16)  and  (17).  To  generalize  (84b)  for  a  system  with  c  servers, 
note  that 


*0  =  p0 

*1  =  pi  +  P2 

n2  ~  P3  +  p4  +  p5 


Vl  ”  Zj=£-C+1 

(85a) 

with 

i  =  c(c+l)/2  -  1  . 

(85b) 

For  c  _<  k  _<  m  , 

_  (c+l)k-i-l+c 
k  Zj=(c+l)k-£-l  pj  ' 

(85c) 
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With  multiple  servers,  the  expected  value  and  variance  of  the  number  of  busy 
servers  (Nbusy)  are: 


and 


^busy1  ”  sk=l  k  "k  +  CZk=c+l  "k 


V[Nbusy]  '  Ek-1  k2"k  +  "k  -  (E[NbUsy)>2  ' 


With  c  servers,  a  queue  forms  when  the  number  of  customers  in  the  system 
N  exceeds  c.  The  number  in  queue 

o  jr  o 


VNsys-c  •  Nsys-C' 


From  this 


m 


N  =  F 
q  Ek=c+1 


(k-c)ir. 


and 


V[Nq]  "  Ek=c+1  (k“c)2irk  “  Nq  • 


(86) 

(87) 

(88) 

(89) 

(90) 


In  all  these  cases  Little's  formula  applies.  Hence,  the  mean  waiting  time 
in  the  system  (W)  and  the  mean  waiting  time  in  the  queue  (W  )  are  obtained 
via  (21)  and  (22),  respectively. 


27.  For  the  specific  case  of  c  equal  to  2  (with  which  we  are  dealing),  the 
conditional  state  probability  density  (q^)  is  obtained  from  p^  via  a 
procedure  similar  to  that  followed  for  the  single-server  case  (pgf.  11). 

In  employing  Bayes1  theorem,  states  having  the  same  number  of  customers 
in  the  system  are  grouped.  Thus, 

q0  -  mp0/K 

qj  =  (m-lJpj/K 
q2  =  (m-l)p2/K  ; 

and,  in  general  for  2  £  k  £  m-1  and  1  £  j  _<  3, 

q3(k-l)+j-l  =  K_1  <m-k>P3(k-l)+M  '  <91> 
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™:  2l°f  Heq“ati0nS.COntrasts  <»>  for  a  Single  server.  The  constant 
K  Obtained  by  summing  over  all  states,  equating  the  sum  to  unity. 

n  the  general  multiple-server  case,  the  conditional  state  p.d.f.  is 
obtained  from 


q0  =  mP[)/K 

qx  =  (m-l)Pl/K 
q2  =  (m-l)p2/K 
q3  =  (m-2)p3/K 
q4  =  (m-2)p4/K 
q5  =  (m-2)p5/K 


...thru  all  states  with  less  than  c  customers  in  the  system.  The 
last  such  state  has  index  £,  where 


or 


Thus, 


4  =  (c-1) (c+2)/2 


*  =  c(c+l )/2  -  1  . 


q„  =  (m-c+1 )p  /K  . 


(92a) 

(92b) 


Additionally,  for  the  states  in  which  the  number  of  customers  (k)  in  the 
system  >  c: 


q j  =  (m-k)pi/K  , 


with 


i  =  (c+1) (k-c)  +  j  +  £  , 

1  1  j  <  c+1  ,  c  <  k  <  m-1  . 

For  2  servers,  the  conditional  probability  that  an 
k  others  in  the  system  is  given  by 

p0  =  q0 
P1  =  ql  +  q2 


(93a) 

(93b) 


arriving  customer  finds 


Pl<  '  >1  q3(k-l)+j-l 


,  2  <  k  <  m-1  . 


(94) 
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This  result  corresponds  to  equation  (26)  for  a  single-server  system.  A 
generalization  of  this  result  to  a  c-server  system  is  the  following: 

p0  =  q0 

P1  =  ^1  +  q2 

p2  =  q3  +  q4  +  q5 

°c-l  ■  Ej=l  Vc+J  •  (95a) 

And, 

_  ^c+1 

pk  j=l  q(c+l ) ( k-c )+£+j  ,  c  _<  k  _<  m-1  .  (95b) 

28.  Three  Servers  With  Erlanq(2)  Service 

Derivation  of  the  state  equations  for  the  three-server  system  with 
Erlang(2)  service  times  follows  the  pattern  used  for  the  two-server  system 
(  PQf •  21  ff).  To  facilitate  recognition  of  states,  double  subscript  notation 
is  also  used  here.  The  state  E..J  signifies  that  i  customers  are  in  the  system, 
with  the  service  stages  for  the  three  channels  occupied  indicated  by  subscript  j. 
The  max  value  of  j  is  c+1  or  4,  in  this  case.  These  states  are  depicted 
iconical ly  and  are  labeled  in  Table  4.  Note  that  the  last  state  for  which 
no  queueing  occurs  is  E^  The  state  transition  diagram  for  this  system  is 
shown  in  Figure  8.  As  indicated  earlier,  one  can  immediately  write  the 
Kolmogorov  equations  by  inspection  from  the  state-transition  diagram. 
Parenthetically,  the  transient  version  of  this  model  was  solved  for  a 

communication  system  having  the  same  transition-state  diagram.  See  p.  14 
of  [5].  The  Kolmogorov  equations,  in  double-subscript  notation,  yield  the 
time-dependent  probabilities  P^ft),  for  admissible  values  of  i  and  j. 

Since  the  steady-state  solution  is  of  interest  here,  time  derivatives  are 
set  to  zero.  Further,  the  development  of  a  matrix-vector  equation  for  the 
steady-state  probabilities  requires  conversion  to  single-subscript  notation. 


[5]  Schlenker,  G.  Op.  Cit.,  May  1977. 
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TABLE  4 


STATES  OF  THE  SERVICE  SYSTEM  WITH 
THREE  SERVERS  AND  TWO  STAGES  OF  SERVICE 


Label 

(i  »j) 

Q 

Configuration 

Stage  1  Stage  2 

Customers  in 

the  System 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1,1 

0 

1 

0 

1 

0 

0 

0 

0 

1,2 

0 

0 

1 

1 

0 

0 

0 

0 

2,1 

0 

1 

0 

2 

i 

0 

0 

0 

2,2 

0 

0 

1 

2 

i 

0 

0 

0 

2,3 

0 

0 

1 

2 

0 

1 

0 

0 

3,1 

0 

1 

0 

3 

1 

0 

1 

0 

3,2 

0 

0 

1 

3 

1 

0 

1 

0 

3,3 

0 

0 

1 

3 

0 

1 

1 

0 

3,4 

0 

0 

1 

3 

0 

1 

0 

1 
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TABLE  4  (Cont) 


Label 

(i  J) 

Q 

Configuration 

Stage  1  Stage  2 

Customers  in 

the  System 

4,1 

1 

1 

0 

4 

1 

0 

1 

0 

4,2 

1 

0 

1 

4 

1 

0 

1 

0 

4,3 

1 

0 

1 

4 

0 

1 

1 

0 

4,4 

1 

0 

1 

4 

0 

1 

0 

1 

For  k  customers  in  the 

system, 

3_<  k  <_  m . 

k,l 

k-3 

1 

0 

k 

1 

0 

1 

0 

k,2 

k-3 

0 

1 

k 

1 

0 

1 

0 

k,3 

k-3 

0 

1 

k 

0 

1 

1 

0 

k,4 

k-3 

0 

1 

k 

0 

1 

0 

1 
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Customers  in 
the  System 
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Figure  8.  State  Transition  Diagram  for  a  Three-Server  System  with  Erlang  (2)  Service  Times 


29.  Steady-State  Equations  (Three  Servers) 

Tie  transformation  from  double  to  single  subscript  notation  in  (73a)  is 
also  used  for  this  case.  Additionally, 

p6  =  P31H 

p7  =  P32^ 

p8  =  P33(a>) 

p9  =  P34^°°^  ’ 

In  general ,  for  3  <_  n  <  m  , 

p4n-6  =  Pn,l(“) 

p4n-5  *  Pn,2^ 

p4n-4  =  pn,3(“) 

p4n-3  =  pn,4(">  '  (96) 

With  these  substitutions,  the  Kolmogorov  equations  in  [5]  become: 

-mXp0  +  yp2  =  0 

mxp0  -  [(m-1) A  +  y]pj  +  yp4  =  0 
yPj  -  [ (m-l)x  +  y]p2  +  2yp5  =  0 
(m-l)Apj  -  [(m-2)A  +  2y]p2  +  yp7  =  0 
(m-l)xp2  +  2yp3  -  [(m-2)A  +  2y]p4  +  2ypg  =  0 
UP4  "  [(m-2)A  +  2y ]pj-  +  3ypg  =  0 
(m-2) Ap3  -  [ ( m- 3 ) A  +  3y]pg  +  yp  =  0 
(m-2)Ap4  +  3yp6  -  [ (rn-3)A  +  3y]p7  +  2yp12  =  0 
(m-2) Ap5  +  2yp7  -  [ (m-3) A  +  3y]pg  +  3ypj3  =  0 
UP8  -  ((m-3)A  +  3y ] Pg  =  0  . 
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For  4  £  n  _<  m-1  , 

(m-n+l)Ap4n_10  -  [(m-n)A  +  3y]p4n_6  +  yp^  =  0 
(m-n+1) xp4n_9  +  3yp4n_6  -  [(m-n)x  +  3y]p4n_5  +  2yp4n  =  0 
(m-n+l)Ap4n_8  +  2yp4n_5  -  [(m-n)A  +  3y]p4n_4  +  3yp4n+1  =  0 
(m-n+l)Ap4n_7  +  yp4n_4  -  [(m-n)A  +  3y]p4n_3  =  0  . 

The  last  four  states  of  the  system  are  described  by: 

Xp4m-10  '  3"P4m-6  “  0 
Xp4m-9  +  3up4m-6  '  3“p4m-5  “  0 
Xp4m-8  +  2pp4m-5  -  3pp4m-4  =  0 

Xp4m-7  +  up4m-4  ’  3up4m-3  =  0  •  <97> 

This  set  of  equations  contains  a  redundancy  since  the  state  probabilities 
must  sum  to  unity.  To  eliminate  this  redundancy  one  can  substitute 

p0  =  1  Ek=l  pk 

into  the  first  equation  in  set  (97),  yielding 
4m  3 

mXZk=l  pk  +  yp2  =  mA  *  (99) 

The  second  equation  in  (97),  which  also  involves  p^,  is  deleted  from  the 
equations  to  be  solved.  The  final  solution  set  has  the  following  matrix- 
vector  form: 

=  b  ,  ( 100a ) 

where  £  and  Jd  are  (4m- 3  X  1)  and  A  is  (4m- 3  X  4m- 3)  . 

In  this  case, 

b'  =  (mA,  0,  0,  ...,  0)  .  (100b) 

These  equations  are  analogous  to  (76),  (77),  and  (78)  in  tiie  two-server 
case.  In  the  general  c-server  case,  the  dimension  of  £  is 

(c+l)m  -  (c-l)c/2  . 
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30.  Statistics  for  the  Three-Server  System 

Recall  that  irk  is  defined  as  the  probability  that  k  customers  are  in  the 
system.  Elements  of  the  system  state  probability  vector  £  determine  the 
elements  of  tt.  The  general  relationship  for  a  system  with  c  servers  is  given 
in  (85).  Particularizing  this  for  3  servers, 

*0  =  po 
*1  =  pi  +  P2 


=  P2  +  P4  +  P5 


and,  for 


*k 


3  _<  k  _<  m  , 

z4k'3  p 
j=4k-6  pk  * 


(101) 


Other  statistics  for  this  case  are  obtained  from  their  c-server  generalizations. 
Specifically,  the  mean  and  variance  of  busy  servers  are  obtained  from  (86)  and 
(87).  The  mean  and  variance  of  number  of  customers  in  the  queue  derives  from 
(89)  and  (90),  with  mean  and  variance  of  number  in  the  system  given  by  (16)  and 
(17).  The  conditional  probability  distribution  for  an  arrival  finding  k  other 
customers  in  the  system  is  obtained  from  (95).  Since  Little's  formula  holds 
for  a  c-server  system,  the  mean  waiting  time  in  the  system  (W)  and  in  the 
queue  (Wq)  are  obtained  from  (21)  and  (22),  respectively.  For  a  system  with 
c  servers,  the  probability  that  an  arriving  customer  must  wait  in  queue  is 

P{arrival  must  wait}  =  z™"* 


or 


1  - 


z 


c-1 

k=0 


Jk  * 


(102) 


Each  of  the  above  statistics  is  calculated  for  c  =  1,  2,  or  3  servers  in 
the  computer  program  FINITE. ME2.Q,  listed  in  Annex  A. 
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ANNEX  A 


COMPUTER  SOURCE  PROGRAMS 


The  computer  programs  in  this  annex  are  used  to  calculate  a  variety  of 
steady-state  statistical  properties  of  queueing  systems  having  a  finite 
customer  population  and  a  multi-server  service  system.  Calculations  are 
made  for  each  of  two  assumptions  regarding  the  probability  distribution  of 
service  times.  The  random  service  time  is  treated  as  derived  from  an  Erlang 
distribution  with  shape  parameter  2  in  the  program  FINITE. ME2.Q.  The  service 
time  is  treated  as  an  exponential  random  variable  in  the  routine  BEST. SERVICE. 
Each  of  these  programs  call  utility  routines  for  inverting  and  multiplying 
matrices:  MAT. INVERSE,  MAT.MPY,  and  MAT.VEC.MPY.  All  of  these  programs  are 
writter  in  SIMSCRIPT  2.5  for  the  PRIME  750  minicomputer.  However,  the  programs 
do  not  employ  features  unique  to  this  computer.  Cross  reference  lists  are 
included  with  program  statements  to  facilitate  the  identification  of  variable 
type  and  to  facilitate  the  location  of  variables  in  a  program. 

A  brief  (50  lines  of  executable  code)  driver  program  TEST. FINITE. ME2.Q  reads 
data  interactively  from  the  terminal  and  echoes  these  inputs.  This  program 
calculates  certain  queueing  statistics  for  a  service  system  with  an  infinite 
population  of  customers  having  Poisson  arrivals  and  with  Erlang(2)  service 
times.  These  quantities  can  be  compared  with  statistics  from  systems  having 
a  finite  customer  population.  The  driver  program  in  turn  calls  FINITE. ME2.Q 
and  BEST. SERVICE.  The  first  of  these  programs  is  largest  with  415  lines  of 
executable  code.  The  statistics  listed  in  the  program's  beginning  comments 
are  optionally  printed  by  the  routine. 

The  routine  BEST. SERVICE  has  174  lines  of  executable  code.  This  program 
has  two  calcul ational  modes:  one  to  produce  queueing  statistics  for  the 
parameters  supplied  and  another  to  sequentially  optimize  the  number  of  servers, 
subject  to  prescribed  constraints.  In  the  present  application  only  the 
calcul ational  mode  is  invoked. 

The  major  problem  segments  of  all  programs  are  announced  via  comment 
statements,  which  are  identified  with  starting  double  quote  marks.  Inputs 
to  the  driver  program  are  read  from  the  terminal  in  response  to  prompting 
messages.  Output  is  sent  directly  to  the  terminal  for  display.  Since  the 
output  is  lengthy,  it  is  recommended  that  a  COMO  file  be  established  to 
display  or  print  it. 
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ANNEX  B 


HETEROGENEOUS  MACHINE  POPULATION 


The  theory  presented  in  the  body  of  this  report  treats  the  case  of  a 
homogeneous  population  of  machines  (customers),  whose  reliability  and 
maintainability  (RAM)  parameters  are  the  same  for  all  machines.  Thus, 
equations  (1),  (2),  and  (3)  do  not  distinguish  between  categories  or  types 
of  machines  having  possibly  different,  intrinsic  mean  times  between  failure 
(MTBF)  and  mean  times  to  repair  (MTTR) .  For  this  ideal  case,  the  approach 
presented  in  the  body  of  the  report  is  exact. 

In  many  instances  the  population  of  machines  is  heterogeneous.  In  this 
case  an  approximation  is  needed  for  the  RAM  parameters  which  are  equivalent 
to  their  homogeneous  counterparts.  This  annex  discusses  an  approximation 
which  can  be  used  in  conjunction  with  the  previous  analytical  results  to 
yield  queueing  statistics  which  are  in  good  agreement  with  simulated  results 
for  tandem  production  operations. 

On  a  production  line  different  machine  operations  may  each  have  several 
machines  of  a  common  type.  However,  each  type  of  machine  (or  operation)  may 
have  RAM  parameters  different  from  parameters  of  other  types.  Furthermore, 
in  this  setting  the  machines  in  a  particular  operation  (n)  do  not  operate 
continuously.  Their  failure  behavior  will  reflect  this  condition.  Assuming 
that  the  machine  operations  are  asynchronous,  a  machine  must  wait:  (a)  if  the 
downstream  buffer  is  full,  (b)  if  the  upstream  buffer  is  empty,  or  (c)  if  the 
machine  is  down  for  repairs.  The  definition  of  arrival  rate  per  machine  (x) 
requires  the  use  of  the  recriprocal  of  the  actual --as  opposed  to  intrinsic-- 
mean  cnock  time  between  failures.  Call  this  the  effective  mean  time  between 
failures  MTBF*.  Similarly,  the  effective  mean  time  to  repair  (MTTR*)  must 
reflect  the  diversity  of  types  of  machines  repaired,  the  MTTR  for  each  type, 
and  the  relative  demand  for  repair  by  type.  Then,  the  repair  parameter  ( y ^ 
used  in  homogeneous  queueing  theory  is  the  reciprocal  of  MTTR*. 

The  following  method  is  a  partial  (and  imperfect)  attempt  to  account  for 
the  non-operating  time,  in  addition  to  maintenance  downtime,  which  the  machines 
of  a  particular  operation  experience.  Generally,  the  production  rates  of 
acceptable  parts  from  the  machine  operations  are  not  identical,  i.e.,  the  line 
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is  unbalanced.  The  machines  of  an  operation  which  does  not  have  minimal 
thruput  must  wait  a  portion  of  the  time.  I  define  the  operational  ratio  of 
operation  n  as  the  fraction  of  the  average  available  machine-time  which 
machines  at  this  operation  must  work  in  order  to  meet  line  thruput.  Reference 
is  made  to  Table  B-l  for  a  complete  list  of  symbols  and  terms  used  in  this 
discussion.  The  average  number  of  machine-minutes  actually  worked  during  a 
day  by  machines  at  operation  n  is,  then, 

U(n)D  A(n)p(n)  ,  (B-l) 

where  the  symbols  U,  D,  A,  and  p  are  defined  in  Table  B-l.  The  expected 
number  of  maintenance  actions  required  at  this  operation  is  this  total 
operating  time  divided  by  the  MTBF(n).  This  number  is  denoted  v(n).  The 
total  of  average  daily  maintenance  actions  required  is  denoted  by  X. 

Formally, 

X  =  zjjsl  v(n)  .  (B-2) 

The  total  daily  average  maintenance  time  required  is  the  following  sum  over 
all  operations: 

T  =  eJJ=1  v(n)MTTR(n)  .  (B-3) 

The  effective  mean  time,  per  machine,  between  failures  is  the  maximum 
daily  machine-minutes  (D  M)  divided  by  X: 

MTBF*  =  D  M/X  .  (B-4) 

The  effective  mean  time  to  repair  is 

MTTR*  =  T/X  .  (B-5) 

The  equivalent  arrival  rate  of  maintenance  actions--l/MTBF*--is  the  sum  of 
required  actions  per  unit  time  over  all  machine  types  divided  by  the  total 
number  of  machines.  The  equivalent  mean  time  to  repair,  or  reciprocal 
equivalent  service  rate,  is  a  weighted  sum  of  MTTR's  for  each  machine  type. 

The  weights  in  this  expression  represent  the  fraction  of  the  total  maintenance 
actions  contributed  by  machines  of  this  (n)  type. 

At  this  point  the  reader  may  have  noticed  that  one  of  the  outputs  of  the 
queueing  model  (W^)  is  needed  to  define  the  equivalent  values  of  X  and  y. 
Availability  (A)  of  machines  of  a  given  type  in  a  manufacturing  operation  is 
used  to  estimate  the  thruput  of  that  operation.  This,  in  turn,  is  used  to 
develop  parameters  x  and  y  of  the  queueing  model.  However,  A  must  consider 
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the  expected  waiting  time  in  the  maintenance  queue  (W^)  as  well  as  the  intrinsic 
RAM  characteristics  of  that  type  of  machine.  Since  A  involves  W  ,  which  is  not 
yet  known,  an  iterative  solution  technique  is  required.  This  process  starts  by 
assuming  that  Wq  is  zero,  calculates  A  for  each  machine  type,  which  leads  to 
estimates  of  X  and  y.  Use  of  X  and  y  in  the  queueing  model,  then,  yields  an 
improved  estimate  of  W^.  The  better  estimate  of  permits  better  estimates 
of  x  and  y,  and  so  forth.  To  accelerate  the  conversion  of  this  iterative  process, 
it  is  helpful  to  use  a  value  of  Wq  (the  same  for  all  machines)  which  is  the 
average  of  the  i  th  and  (i-1)  th  iterates.  Convergence  in  some  instances  is  not 
uniform.  Computational  experience  indicated  that  truncation  of  this  iterative 
process  occurs  rapidly  using  the  following  convergence  criterion:  the  relative 
difference  in  values  of  for  the  last  two  iterations  must  be  less  than  2%. 

Sound  computer  coding  practice  provides  a  sure  escape  from  the  iteration  loop 
after  a  maximum  number  of  iterations  has  occurred.  In  the  TANDEMT  program, 
the  maximum  was  set  to  20,  altho  iterative  truncation  has  yet  to  occur. 
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TABLE  B-l 


> 


NOTATION  USED  IN  OBTAINING 

EFFECTIVE  ARRIVAL  AND  SERVICE  RATES 

Symbol 

Description 

D 

length  of  workday  (minutes) 

U(n) 

the  number  of  (units  of)  machines  of  a  common  type  in 
operation  n 

N 

the  number  of  machine  types  or  machine  operations 

M 

N 

the  total  number  of  machine  units  ,  =  En  U(n) 

MTBF(n) 

the  mean  operating  time  between  failures  (minutes) 
of  machines  of  type  n 

MTTR(n) 

the  mean  time  (minutes)  to  repair  a  machine  of  type  n, 
given  maintenance  resources 

WQ 

the  mean  time  (minutes)  a  machine  spends  in  the  repair  queue 

A  (n ) 

the  steady-state  availability  of  machines  of  type  n, 

=  MTBF(n)/[MTBF(n)+MTTR(n)+WQ] 

R(n) 

unit  machine  rate  (per  min)  at  operation  n 

AP(n) 

average  daily  production  capacity  at  operation  n 
=  U(n)D  R(n)A(n) 

a(n) 

acceptance  rate  of  parts  from  machines  at  operation  n  only 

C(n) 

r 

K 

cumulative  acceptance  rate  of  parts  passing  operation  n 
per  part  starting  the  first  operation 

*  "3=i  •«> 

CP(n) 

average  daily  production  capacity  of  acceptable  parts 
at  operation  n,  =  AP(n)  a(n) 

n* 

value  of  n  for  which  CP(n)  is  a  minimum 

CP(n*)  <_  CP(n)  ,  1  <  n  <N 

RP(n) 

average  daily  parts  production  at  operation  n  required 
to  meet  the  average  daily  line  thruput,  =  AP(n*)  C(n)/C(n*) 
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TABLE  B-l  (Cont'd) 


NOTATION  USED  IN  OBTAINING 
EFFECTIVE  ARRIVAL  AND  SERVICE  RATES 


Symbol 

p(n) 

v(n) 

X 

T 

MTBF* 

MTTR* 

X* 


Description 

operational  ratio  for  machine  operation  n,  =  RP(n)/AP(n) 

average  daily  number  of  maintenance  actions  for  machines 
of  common  type  in  operation  n,  =  U(n)D  A(n)p(n)/MTBF(n) 

total  average  daily  maintenance  actions  required, 

=  E^  v(n) 
n  '  ' 

total  average  daily  maintenance  time  (min)  required, 

=  zJJ  v(n)MTTR(n) 

effective  mean  time  (min)  between  failures,  =  D  M/X 
effective  mean  time  (min)  to  repair  all  machines,  =  T/X 
effective  arrival  rate  (1/minute),  =  1/MTBF* 
effective  service  rate  (1/minute),  =  1/MTTR* 


* 
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